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ABSTRACT 

We present a new iterative method for constructing equilibrium phase models of stel- 
lar systems. Importantly, this method can provide phase models with arbitrary mass 
distributions. The method is based on the following principle. Our task is to gener- 
ate an equilibrium iV-body system with a given mass distribution. For this purpose, 
we let the system reach equilibrium through its dynamical evolution. During this 
evolution we hold mass distribution in this system. This principle is realized in our 
method by means of an iterative procedure. We have used our method to construct 
a phase model of the disc of our Galaxy. In our method, we use the mass distribu- 
tion in th e Galaxy as input data. Here we used two Galactic density models (sug- 
gested bv lFlvnn. Sommer-Larsen fc Christensen|[l996HDehnen fc Binnev|[l99 8a). For 
a fixed-mass model of the Galaxy we can construct a one-parameter family of equi- 
librium models of the Galactic disc. We can, however, choose a unique model using 
local kinematic parameters that are known from Hipparcos data. We show that the 
phase models constructed using our method are close to equilibrium. The problem of 
uniqueness for our models is discussed, and we discuss some further applications of 
our method. 
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1 INTRODUCTION 

The construction of equilibrium phase models of galaxies is 
an important area of research in galactic astronomy. Such 
models are of interest from a number of points of view. They 
are important for understanding the dynamics of galaxies, 
and they are necessary for defining the initial conditions in 
iV-body models of stellar systems. 

In this paper, we consider the problem of constructing 
an equilibrium model of the stellar disc of a spiral galaxy. 
Our purpose is to construct a phase model of the Galac- 
tic stellar disc. Various approaches to solving this prob- 
lem have been suggested (s ee, for example, the review in 
iRodionov fc Sotnikovalkood hereafter RS06). 

The first approach is based on Jeans equations (equa- 
tions for moments of the equilibrium velocity distribution 
function). One such met hod of constructin g equilibrium disc 
models was described bv lHernquistl (Il993l ). An advantage of 
this method is its relative simplicity. It is applicable for a 
stellar disc with an arbitrary density profile and any external 
potential. It has, however, a significant drawback. The sys- 
tem of Jeans equations used is not closed, so it is necessary 
to introduce an additional condition in order to close it. As 
a result, the constructed model is often far from equilibrium, 
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as we showed when we used the closure condition suggested 
bv lHernquistl l| 19931) . A more detailed critical analysis of this 
method is given in RS06. 

The second approach is based on Jeans theorem, accord- 
ing to which any function of motion integrals is a solution 
of the st ationary collisionless Boltz mann equation (see, for 
example. iBinnev fc Tremaine|[l987h : that is, it is an equi- 
librium distribution function. There is, however, one signif- 
icant disadvantage of such an approach to constructing a 
three-dimensional equilibrium model of a stellar disc. Two 
integrals of motion are well known for axisymmetric mod- 
els: E is energy and L z is the angular momentum about the 
symmetry axis. However, for systems having phase density 
f(E, L z ), the dispersions of the residual velocities in the 
radial and vertical directions have to be the same, which 
is in disagreement with observations of spiral galaxies, in 
particular for the solar neighbourhood (see, for example, 
iDehnen fc Binnevll9 98b). Axisymmetric models with differ- 
ent velocity dispersions in the radial and vertical directions 
may be constructed if the phase density depends on three 
integrals of motion f(E, L z , I3), where 73 is the third inte- 
gral of motion. However, an expression for the third integral 
is not known for the general case. It is possible to use the 
energy in vertical oscillations as the third integral when the 
residual velocities are much lower than the rotation velocity 
with respect to the symmetry axis (cold thin disc) . In such a 
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way, one can construct models of appro ximately exponential 
stellar discs, as done, for example , by iKuiiken fc Dubinskil 
(1995): lwidrow fc DubinskH (2005). These authors also de- 
scribe the procedure of phase density construction for mul- 
ticomponent models of disc galaxies. 

One further original metho d for construct i ng p hase 
galactic models was developed by ISchwarzschiidl (|l979h . In 
this method, it is assumed that the total galactic potential 
is known. A large number of orbits (library of orbits) in 
this potential are constructed, and a model consisting of the 
particles placed on these orbits is constructed in such a way 
that the resulting model has an initial density profile. We 
note that this approach is similar to our Orbit. NB method 
described below. 

In this paper, we use a new iterative method proposed 
in RS06. In RS06, the iterative method was applied to con- 
struct a model of the stellar disc of a spiral galaxy (the 
problem under consideration). In our notation, this real- 
ization of the iterative method is termed Nbody.SCH. It 
has, however, a number of disadvantages. The Nbody.SCH 
method has a problem with the construction of a relatively 
cold stellar disc: using the Nbody.SCH method it is not 
possible to construct a sufficiently cold equilibrium stellar 
disc. The stellar disc of our Galaxy is rather cold. The 
model of the stellar disc of our Galaxy constructed using the 
Nbody.SCH method is notably far from equilibrium. More- 
over, the Nbody.SCH method cannot be directly applied to 
a stellar system with arbitrary geometry (for example, it 
cannot be applied to elliptical galaxies). 

The first objective of our work is to develop a new real- 
ization of the iterative method without the disadvantages of 
the Nbody.SCH approach. The second objective is to con- 
struct a phase model of the Galactic disc using this new 
method. 

Using the iterative method we can construct a phase 
model of a stellar system with a given mass distribu- 
tion. Therefore, in order to construct a phase model 
of the Galactic disc we need the mass distribution of 
the Galaxy. A number of density models for the Milky 
Way Galaxy have been constructed. We use only two 
of them (see Flvnn. Sommer-Larsen fc Christensen] 1 19961 ; 
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The density models used are described in Section [2] In 
Section [3] we present the iterative method and its modifica- 
tions. Two versions of the iterative method are considered in 
detail. In our notation, these methods are called Orbit. NB 
and Nbody.NB. The Orbit. NB method gives fairly specific 
and probably non-physical models. Although such models 
are probably non-physical, the fact that such equilibrium 
models exist is of interest. It gives some insight into the 
problem of uniqueness of phase models of the stellar disc. 
The models constructed by the Orbit. NB method are dis- 
cussed in Section [4] We think that the models constructed 
by Nbody.NB method are physical. In Section we present 
phase models for the Galactic disc constructed using this 
method. It is important, that Nbody.NB solves the first ob- 
jective of our work (see above). A summary of the results is 
given in Section [6] wherein we also discuss further applica- 
tions of our method and models. 



Galactic density models in 
have chosen two of them 



There are many 
the lite rature. We 

(see Flvnn. Sommer-Larsen fc Christensen! Il99i 

IDehnen fc Binned 1 1998al ). We note that IDehnen fc Binne 
(1998a) presented a whole family of density models. We 
have chosen the model '2' from this paper. Both models 
under consideration are axisymmetric. 

Let us briefly outline the models used. 



2.1 Model of 

iFlvnn. Sommer-Larsen &: Christensenl (|l99f 



This model contains three main components: a dark halo, 
central component, and disc. For the dark halo, the authors 
used a logarithmic potential fsee lBinnev fc Tremainelll987l . 
p. 46): 



(1) 



where Vh and rg are the halo parameters (circular velocity 
at large r and halo scalelength) , R is the cylindrical radius, 
r — \J R? + z 2 is the spherical radius. 

The central component consists of two spherical sub- 
systems. The first one represents the bulge+stellar halo; 
the second one is the inner core of the Galaxy. Each 
component is approxima ted by a Plummer sphere (see 
iBinnev fc Tremaindll987l . p. 42-43). The expression for the 
whole potential of the central component has the form 



$c(R,z) = - 



GM Cl 



GMc 2 



\Jr 2 + rd 2 y / r 2 + rc 2 



(2) 



where G is the gravitational constant; Mc 1 and rc 1 are the 
mass and scalelength for the first subsystem; and Mc 2 and 
rc 2 are the same parameters for the second. 

The disc in th i s mod el is the superposition of three 
iMivamoto fc Nagail l|l975h discs. The whole disc potential 
has the form 



<f D (R,z) = -J2 



GMd„ 



T[ y/R 2 + (a n + Vz 2 + b 2 ) 2 



(3) 



Here, b is the disc scaleheight (which is the same for all three 
components); the parameters a n are the disc scalelengths; 
and the values Mp. are the masses of the disc co mponents. 

IFlvnn. Sommer-Larsen fc Christensenl l|l996h have sug- 
gested the following values for the above parameters: 
r H = 8.5 kpc, Vh = 210 kms -1 , 
r Cl = 2.7 kpc, M Cl = 3.0 • 10 9 M a , 



rc 2 = 0.42 kpc, M c , 

b = 0.3 kpc, 

M Dl = 6.6 ■ 10 10 M Q 



1.6 ■ 10 10 Mi 



oi = 5.81 kpc, 



) 10 

Md 3 = 3.3 • 10 9 M Q , a 3 = 34.86 kpc. 

Table 1 in IFlvnn. Sommer-Larsen fc C hristensen (1996) 
contains a small misprint: instead of Vh = 220 kms -1 it 
should be Vh = 210 kms -1 (Flynn, private communica- 
tion). We note that one of the disc components (n = 2) 
has a negative density; however, the total density in the 
disc is positive, and the model is therefore physical. Over a 
large range in R, the disc density profile is approximately 
exponential, with a scalelength of approximately 4 kpc. This 
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Table 1. 


Parameters for 


three disk 


components 


in the 


DB2 model. 










comp. 


E d , Mope" 2 


R d , kpc 


Rm, kpc z d , kpc 


thin disk 


1022 


2.4 





0.18 


thick disk 


73.03 


2.4 





1 


ISM 


113.6 


4.8 


1 


0.04 



Table 2. Parameters for the two spheroidal components in the 
DB2 model. 



comp. 


P0, M pc 3 


ro, kpc 7 


(3 


H 


ft, kpc 


bulge 


0.7561 


1 1.8 


1.8 


0.6 


1.9 


halo 


1.263 


1.090 -2 


2.207 


0.8 


CO 



is possibly an overestimated value (Flynn, private commu- 
nication). Hereafter, we will refer to this model as FSLC. 
Fig. [T] shows the dependences of cumulative masses M(r) 
and circular velocity curves for the whole FSLC model, for 
all components without the disc, and for the disc only. 

2.2 Model of lDehnen &: Binnevl (|l998ah 

In addition to the FSLC model, we consider one m odel of 
the family suggested by iDehnen fc Binnevl ljl998al ). Every 
model of this family consists of five components. There are 
three disc components (interstellar medium (ISM), thin stel- 
lar disc, and thick stellar disc) and two spheroidal compo- 
nents (dark halo and bulge). 

The distribution of volume density in each disc compo- 
nent has the form 



pd{R,z) = 7^- exp 



Rm 
Rd 



R_ 

Rd 



\z\ 

Zd 



(4) 



Here T, d is the central surface density of the component, and 
parameters R d and Zd give scalelength and scaleheight of the 
component. By using the parameter R m one can introduce 
a central density depression in the disk. 

The density distribution for each spheroidal component 
has the form 



ps(m) = po 



where 



ro 



m 
ro 



7-/3 



^R? + q~ 2 z 2 . 



exp(-m 2 /V t 2 ), (5) 



(6) 



Here po, ro, 7, (3, q, r t are the parameters of the spheroidal 
components. 

We use model 2 from this paper, hereafter the DB2 
model. This choice is somewhat arbitrary. We do not con- 
sider other models from this paper because a comparison 
of the different Galactic models is outwith the goals of this 
paper. The parameters of the DB2 model are shown in Ta- 
bles Q] and El The details of the const ruction of this model 



are given in Dehnen fc Binnevl (l99 8a). 



Here we construct an equilibrium TV-body model of the 



stellar Galactic disc. As the stellar disc, we take only the 
thin stellar disc from the DB2 model. The construction of a 
two-component stellar disc in this model will be the subject 
of future investigations. 

In addition to the density, we need to calculate 
the potentials of the various components in the DB2 
model. The potentials were numerically calculated using 
the code GALPOT by Walter Dehne n The method of 
potent ial determination is described in (|Dehnen fc Binnevl 
1998a). The code was take n from the NEMO package 
( |http:/ /astro. udm.edu/nemo| lTeubenlll995l ). 

The cumulative mass profile M(r) and circular velocity 
curve for the DB2 model are shown in Fig. [T] for the whole 
model, for all components apart from the thin stellar disc, 
and for the thin stellar disc only. 

Fig. Q] shows that the FSLC and DB2 models are differ- 
ent. The FSLC model has a more massive and concentrated 
bulge with respect to the DB2 model. In particular, this 
massive bulge is the reason for the central peak in the ro- 
tation curve for the FSLC model. Furthermore, the relative 
disc contribution in the whole mass and the circular velocity 
curve for the inner model (R < 8 kpc) is significantly higher 
in the DB2 model than in the FSLC model. 



3 ITERATIVE METHOD FOR EQUILIBRIUM 
MODELS CONSTRUCTING 

3.1 Basic idea of iterative method 

The iterative method is used to construct TV-body models 
close to equilibrium and with a given mass distribution (see 
RS06 for details). The basic idea of this approach is as fol- 
lows. First, we generate an TV-body system with a given 
mass distribution but with arbitrary initial particle veloci- 
ties (which, for example, can be taken as zero). Furthermore, 
we let the system reach equilibrium through its dynamical 
evolution. During this evolution we hold mass distribution 
in the system. If necessary, some parameters of the velocity 
distribution can be fixed. This is achieved in the following 
way. 

The general algorithm of the iterative method is as fol- 
lows. 

(1) An N-body system with a given mass distribution but 
with arbitrary particle velocities is constructed. The veloci- 
ties can, for example, be taken to be zero. 

(2) The system is evolved on a short time-scale. 

(3) We construct a new N-body system, with the same 
given mass distribution but with velocities chosen according 
to the velocity distribution in the system already evolved. 
We note that, if there are some limitations on the veloc- 
ity distribution, this distribution should be corrected taking 
into account these restrictions (see the discussion below). 

(4) We return to point 2. We repeat such cycles until the 
velocity distribution stops changing. 

As a result, we obtain an TV-body model close to equi- 
librium that has a given density profile (see RS06 and our 
results below for details). 

We can discuss the iterative approach in a more general 
manner. When one needs to find an equilibrium state of an 
arbitrary dynamical system, but so that this state has some 
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Figure 1. The dependence of cumulative masses on radius (left panels) and the circular velocity curves (right panels) for various 
components in the density models FSLC and DB2. Here the cumulative mass M(r) is the mass inside the sphere of radius r. Solid lines 
show the dependences for a whole model; long-dashed lines correspond to the whole model without the disc; and short-dashed lines 
correspond to the stellar disc only. We use the thin stellar disc in the DB2 model as the stellar disc. 



necessary properties (in the case under study, the dynami- 
cal system is a set of gravitating points and the necessary 
property is the density profile), one can simply give a possi- 
bility for the system itself to tune to the equilibrium state, 
holding the necessary parameters. 

The idea of our iterative method is simple. Its realiza- 
tion in practice is more complicated. The main difficulty is 
the third stage, when it is necessary to construct a model 
with the same velocity distribution as the evolved model 
from the previous iteration step. 

Below we discuss an application of the iterative proce- 
dure to the problem of constructing an equilibrium model 
of the Galactic disc. 



3.2 Equilibrium models of stellar disks 

3.2.1 Family of equilibrium models 

Our task is to construct an equilibrium model of the stellar 
disc of our Galaxy. We consider all Galactic components as 
axisymmetric, and can formulate our task in the following 
way. We need to construct an equilibrium iV-body model of a 
stellar disc with a fixed density distribution pd\ B k(R,z) that 
is embedded in the rigid external potential & C xt(R, z), where 
$ext {R, z) is created by all Galactic components except the 
stellar disc (i.e. the dark halo and bulge). 

It can be expected that at least a one-parameter fam- 



ily of equilibrium models will exist when the functions 
Pdisk(-R, z) and $ext(-R, z) are fixed. The parameter of this 
family is the fraction of kinetic energy contained in the resid- 
ual motions. 

The reason why this family exists is as follows. It is 
possible to show that, if p<jisk(^, z) and <3> e xt(fi, z) are fixed, 
then for all equilibrium discs the total kinetic energy should 
be the same. This is a direct consequence of the virial the- 
orem. This kinetic energy can, however, be distributed be- 
tween regular rotation and residual velocities in different 
manners. Cold equilibrium models exist when a large frac- 
tion of the kinetic energy is concentrated in the regular ro- 
tation (an extreme case is the model with circular orbits), 
and hot equilibrium models exist when a large fraction of 
the kinetic energy is concentrated in the residual motions. 

In RS06, the authors showed that, for exponential disc, 
a one-parameter family of models can be constructed by the 
iterative method for fixed pdisk(-R, z) and <J> cx t(-R, z). If the 
iterations are started from different initial states, different 
models result; however, they all form a one-parameter fam- 
ily. As expected, the parameter is the fraction of kinetic en- 
ergy concentrated in the residual motions. In order to obtain 
a definite model from this one-parameter family, the method 
suggested in RS06 can be used. We simply fix the fraction 
of kinetic energy during the iterative process. In principle, 
we could fix any value characterizing the "heat" degree of 
the disc. The authors of RS06 suggest using the value of an- 
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gular momentum about the z-axis (symmetry axis) as this 
parameter: 



,Ri ■ 



(7) 



where m,, v v i, Ri are the mass, azimuthal velocity, and 
cylindrical radius of the i-th particle. 

We fix the value of L z at each iteration step. When we 
have constructed a new model (which has the same velocity 
distribution as the slightly evolved model from the previous 
iteration step), we correct the azimuthal velocities so that 
the total angular momentum of the system is the same as 
the fixed value of L z . This is done as follows. Let L z be this 
fixed angular momentum, and let L' z be the current value 
of the angular momentum. The new azimuthal velocities of 
particles are prescribed as follows: 



(L z -L' z ) 
CL R i m i 



(8) 



where v' vi is the current value of the azimuthal velocity of 
the i-th particle, and v v i is the corrected azimuthal velocity 
of the i-th particle. 

We note that by using the iterative method with fixed 
L z , one can construct colder models with respect to the 
ones without L z fixed (because the cold stellar disc tends 
to heating during the dynamical evolution) . The stellar disc 
of the Galaxy is rather cold. Thus it is difficult to construct 
a cold model of the Galactic disc without fixing L z . We 
therefore fix L z in all our models. 



3.2.2 Variants of the velocity distribution "transfer" 

We first discuss the core of the iterative method, namely an 
algorithm to transfer the velocity distribution (item 3 in the 
iterative procedure). 

The transfer problem is as follows. We have an "old" 
model. This is an evolved model from the previous iteration 
step from which we would like to copy a velocity distribu- 
tion. Moreover, we have a "new" model that is constructed 
according to the fixed density distribution. We have to give 
the velocities to the particles in the new model using the 
velocity distribution in the old model. How do we do this? 

In RS06, the authors used an algorithm of velocity dis- 
tribution "transfer", which is based on assuming that the 
particles have a truncated Schwarzschild velocity distribu- 
tion. We describe this approach briefly. We take a disc model 
(old model) from which we are going to "copy" the veloc- 
ity distribution. The model is divided along the _R-axis into 
various regions (concentric cylindrical tubes). For each re- 
gion, we calculate four velocity distribution moments (v v , 
o~R, cr^, a z — the mean azimuthal velocity and three dis- 
persions of residual velocities along the directions R, ip, and 
z). These moments are used for the velocity choice in the 
new model. We assume that the velocity distribution is the 
Schwarzschild one, but without the particles that can go out 
of the disc (see RS06 for details) . 

We slightly modified this scheme of velocity transfer. 
The model is divided into regions not only along the f?-axis, 
but also along the z-axis. The regions are chosen in such a 
way that they all contain similar numbers of particles. 



This method of velocity distribution transfer has, how- 
ever, two drawbacks (even in its modified form). First, an a 
priori assumption is made that the velocity distribution is 
the Schwarzschild one. Second (and more importantly), this 
method cannot be used for other geometry systems (e.g. tri- 
axial elliptical galaxies). 

We have developed another method of velocity distri- 
bution transfer. We believe that it is a more general and 
simpler method. The basic idea of this new method is as fol- 
lows. We prescribe to the new-model particles the velocities 
of those particles from the old model that are nearest to the 
ones in the new model. 

The simplest (although not quite successful, as we show 
below) implementation of this idea is evident. One can pre- 
scribe to each particle in the new model the velocity of the 
nearest particle from the old model. Let us formulate this 
proposition more strictly. For each ith particle from the new 
model, one finds the old-model jith particle with the mini- 



mum value of |r" 



Here, r" ew is the radius vector 



of the ith particle in the new model, and is the radius 
vector of the jih particle in the old model. One then takes as 
the velocity of the ith particle in the new model the velocity 
of the jr'th particle in the old model. 

This simple algorithm has, however, one essential prob- 
lem. If the numbers of particles in the old and new models 
are the same then only about one-half of the particles in the 
old model participate in the velocity transfer. The reason is 
that many old-model particles have a few particles in the 
new model to which they transfer their velocities. At the 
same time, almost one-half of the particles in the old model 
do not transfer their velocities at all. This means that a sig- 
nificant amount of information on the velocity distribution 
will be lost in the transfer process. The noise will therefore 
grow, and this is indeed observed in numerical simulations. 
Thus, one cannot construct an iV-body model close to equi- 
librium by the iterative method described above using this 
transfer algorithm. 

It is, however, possible to modify the transfer scheme in 
order to overcome this failure. An input parameter of this 
improved algorithm is the "number of neighbours" n nb . For 
each particle in the old model, we introduce the parame- 
ter n uae , which denotes the number of times this particle is 
used for velocity copying. At the beginning of the transfer 
procedure, we assume n use = for each particle in the old 
model. We then consider all particles in the new model. For 
each particle from the new model, we find the nearest n n j, 
neighbours in the old model (in this case, the closeness is un- 
derstood as the minimum distance between the particles in 
the old model and the point at which the particle of interest 
from the new model is placed) . We then reveal a subgroup of 
particles that have a minimum n U3e among these n„b neigh- 
bours, and among this subgroup we find the particle that 
is the closest one to the position of the new-model particle. 
We prescribe to the new-model particle the velocity of the 
found particle from the old model. Moreover, we add one to 
the parameter n uae of this old-model particle. 

We note that this algorithm is the same as the previous 
one if n n j, = 1. As mentioned above, the problem with this 
algorithm is that only one-half of the particles take part in 
the velocity distribution transfer. If we take n nb = 10, only 
a small fraction (a few per cent) of old-model particles do 
not take part in the transfer process. Using this improved 
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transfer method in the iterative procedure gives good re- 
sults in the sense that the constructed models are close to 
equilibrium. 

In this method, one can take into account that the 
galactic models under study are axisymmetric. One just has 
to redefine the distance between the particles. That is, one 
can search for the nearest particles in two-dimensional space 
Rz instead of three-dimensional space XYZ. In this case, 
one should transfer the velocities in cylindrical coordinates 
in order to remove any dependence on the azimuthal angle. 
This guarantees that the constructed model has an axisym- 
metric velocity distribution. We use this method when con- 
structing the model of the Galactic stellar disc (see below). 

We note that this method of velocity transfer is univer- 
sal, and can be applied in systems with arbitrary geometry. 
The models constructed in this way are close to equilib- 
rium, but, partly because of this, the method has a small 
disadvantage. The iterations converge much more slowly 
than the ones in the above-mentioned method based on the 
Schwarzschild velocity distribution. The reason for this slow 
convergence is that even intermediate models are fairly close 
to equilibrium, so that the models change only slightly over 
one iteration, and the iterations converge slowly. 

3.2.3 Different ways of system evolution simulations 

We discuss one more way to modify the iterative method. In 
the general algorithm of the iterative procedure, item 2 al- 
lows the model to evolve over a short time-scale. This means 
that one needs to simulate the self-consistent TV-body evo- 
lution of the stellar disc during a short time in the field 
of the external potential "Joxt- There is, however, another 
possibility. Instead of simulating the self-consistent dynam- 
ical evolution of the TV-body system, one can simulate the 
motions of TV massless test particles in the regular galac- 
tic potential $disk + $cxt, where $disk is the disc potential 
corresponding to the density pdisk- We note that the simu- 
lation of a system of TV test particles in a rigid potential is 
much less cumbersome than simulation of the self-consistent 
TV-body system. 

At first glance, both methods have to give practically 
the same results, because the initial stellar disc has the den- 
sity profile pdisk that creates the potential <E>disk. It is ex- 
pected that the iterations will converge to an equilibrium 
state. Therefore in the self-consistent case, the disc at the 
late stages of the iterative process will be close to equilibrium 
and will not significantly change its density profile during a 
single iteration (especially because we follow the evolution 
over a short time-scale). 

However, the iterative methods using these two modes 
of calculation may lead to essentially different results (see 
below). 

3.2.4 Comparison of different realizations of iterative 
method 

We thus have two ways to perform the velocity distribution 
transfer. The first one is based on calculations of the velocity 
distribution moments and assumes a Schwarzschild velocity 
distribution (hereafter we refer to this method as SCH). The 
second is based on prescribing to the particles in the new 



model the velocities of the nearest particles in the old model 
(hereafter we refer to this method as NB). 

We also have two ways of performing system simula- 
tions. The first one is the self-consistent simulation of the 
N-body gravitating system evolution (hereafter we refer to 
this approach as Nbody). The second is the calculation of 
massless particle orbits in the regular potential "l?disk + $oxt 
(hereafter we refer to this approach as Orbit). 

We thus have four versions of the iterative method: 
Nbody.SCH, Nbody.NB, Orbit.SCH, and Orbit.NB. Our 
task is to choose from among them the method that will 
be used to construct a phase model of the Galactic stellar 
disc. 

We will show in Section [4] that the Orbit.NB approach 
gives models that are probably non-physical. However, the 
very fact of the existence of such "strange" equilibrium mod- 
els is of interest, and gives food for thought (see Section [4] 
for details). 

We thus need to choose from among three approaches: 
Nbody.NB, Nbody.SCH, and Orbit.SCH. Our test simula- 
tions have shown the following. The models of hot discs 
constructed by all these approaches are similar. The only 
exception is that the models constructed by the Nbody.NB 
method are slightly closer to equilibrium. However, the mod- 
els of cold discs are significantly different. In particular, this 
concerns models of the Galactic disc, because the Galactic 
disc is rather cold. The models of cold discs constructed with 
the Nbody.NB approach are close to equilibrium, whereas 
the ones constructed with the Orbit.SCH and Nbody.SCH 
approaches are quite far from equilibrium. We can therefore 
conclude that the Orbit.SCH and Nbody.SCH approaches 
are not suitable for constructing equilibrium models of cold 
stellar discs. Moreover, from a methodical point of view, it 
is more correct to use the NB transfer approach, because 
here we do not make any a priori assumptions concerning 
the form of the velocity distribution (in contrast to the SCH 
method) . 

We shall thus use the Nbody.NB approach to construct 
phase models of the Galactic stellar disc (see Section [SJ. 

3.2.5 Technical comments 

We note a few important technical details. 

In the iterative method, there is one parameter — the 
time interval ti of each iteration. This is the time interval 
over which the system evolves during each iteration. How 
do we choose the value of ti? This time should not be too 
short, because in this case the system would have no time 
to evolve during one iteration step. On the other hand, it 
should not be too long. At least, this time should be much 
shorter than the typical times of development of various in- 
stabilities; otherwise, these instabilities may change the sys- 
tem substantially. We cannot suggest any strict criterion for 
choosing ti. We have chosen a value from numerical simula- 
tions. Our simulations have shown that, if we take ti within 
reasonable limits (not too short and not too long), the con- 
structed models are the same (within the noise limits). Of 
course, this is valid only when we use the model with the 
same L z (see Section f3. 2. In any case, a basic test of ev- 
ery method to construct the equilibrium models would be a 
numerical check that the model is close to equilibrium. 

We have used the following modification of the iterative 
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method in several simulations. We do not take a fixed iter- 
ation time, but choose this time randomly within the range 
(0, t™ ax ). If one makes the iterations with the fixed U, the 
following situation may occur in principle. The iterations 
could converge to an artificial non-equilibrium state when 
the model has large changes at intermediate times within 
one iteration, however in the end of iteration, it could have 
the same state as in the beginning of the iteration. Another 
situation that could occur is that the model could jump from 
one state to another, i.e. oscillations between two states oc- 
cur. The random choice of the iteration step enables us to 
avoid such situations. However, if we consider the Nbody.NB 
approach, fixed and random iteration steps give practically 
the same models as output. 

In many of our simulations, we used the following 
method to decrease the CPU time. Initially, we make the it- 
erations have a low accuracy (for example by using a smaller 
number TV of bodies) and then gradually increase the accu- 
racy up to a pre-defined limit. This scheme was used in all 
simulations of Section [5] 

In all our TV-body simulations ( self-consistent 
schem e) , we used the TREE code |Barnes fc Hutl 
1986) and a few other codes fro m the NEMO package 
(http://www.astro.umd.edu/nemo; Tcubcn 1995). We used 



our original codes to simulate the motions of particles in 
the rigid potential. 



4 NON-PHYSICAL MODELS. HYPOTHESIS 
OF UNIQUENESS 

4.1 Models constructed with the Orbit. NB 
approach 

In this section we consider the models constructed with the 
Orbit. NB approach. We will show that these models are 
probably non-physical models. However, the very fact of the 
existence of such models is of interest. It gives some insight 
into the problem of uniqueness of phase models of the stellar 
disc. 

One noteworthy feature of the iterative Orbit. NB ap- 
proach is that usually at fixed pdisk(-R, z), 4> ex t(-R, z), and 
L z essentially different models can be constructed by means 
of this approach, although all other versions of the iterative 
method give similar (within the limits of noise) final models 
at a fixed value of L z . Another important feature of this ap- 
proach is that the velocity distributions of the final models 
are very different from the Schwarzschild distribution. We 
consider this fact in more detail below. 

Let us consider a model constructed by the Orbit. NB 
approach. We take the FSLC density model (see Sec- 
tion I2.1|l ; that is, the disc density pdisk is taken from the 
FSLC model, and the rigid potential <E> ext is generated by 
all FSLC model components except the disc. The disc den- 
sity is taken to be zero at R > i? ma x = 30 kpc or \z\ > 
Zmax = 10 kpc. We take the cold initial model in which all 
particles move along circular orbits. We used 1000 iterations 
for TV = 200, 000, and then 100 iterations for TV = 500, 000. 
The integration step was taken as dt — 0.5 Myr. The num- 
ber of neighbours for the velocity distribution transfer was 
taken to be n n b = 100. We used a scheme with randomized 
iteration time (see Section [3~23]l with if 13 * = 100 Myr. Dur- 
ing the iterative process, we fixed the angular momentum as 




Figure 2. Dependences of v v , ctr, a v , and cr z on R in the 
FSLC.O model constructed by the Orbit. NB method. 



L z = 6.412- 10 13 M Q kpckms" 1 . We note that in all simula- 
tions the following system of units was used: gravity constant 
G = 1, length unit u r = 1 kpc, time unit ut — 1 Myr. In 
this system of units, the chosen value is L z — 0.295 (here- 
after we indicate the parameter L z in this system of units). 
We refer to this model as FSLC.O. We consider only this 
model; however, we emphasize that very different models 
may be constructed by the Orbit. NB approach at fixed L z 
if we take different initial models. 

In addition to the constructed FSLC.O model, we con- 
sider its changes during further dynamical evolution (self- 
consistent TV-body evolution of the constructed stellar disc 
in the field of the external potential <E> ex t). The parame- 
ters of simulation are taken as follows: number of bodies 
TV = 100 000, integration time step dt = 0.04 Myr, softening 
length e = 0.025 kpc. Two last para meters were chosen ac- 
cordin g to the recommendations of iRodionov fc Sotnikoval 
(|2005h . 

The radial profiles for four moments of the velocity dis- 
tribution v v , (tr, a v , and a z , are shown in Fig. [2] It can be 
seen that the profiles of v v , or, and a v have unusual forms. 
They have various peaks and dips. It would seem that such a 
complicated system cannot be stable; however, the FSLC.O 
model is close to equilibrium! When we follow its evolution, 
it emerges that the constructed model conserves structural 
and dynamical parameters very well (see Fig. [3}. 

An interesting question arises in connection with the 
FSLC.O model equilibrium: how do the moments of the ve- 
locit y distribution satisfy the equilibrium Jeans equations 



(see iBinnev fc Tremain 
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Here $tot = 5>disk + 3>oxt- Fig. [4] shows the radial profiles of 
v v , o v , and cr z from the FSLC.O model and from the Jeans 
equations ((9]) (see also RS06) . It can be seen that the model 
follows the Jeans equations very well. This is an unexpected 
finding, especially considering such unusual moment profiles. 
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Figure 3. Initial evolutionary stages for the FSLC.O model. The upper snapshots show the disc views face-on for various moments 
of time (0, 100, 200, 300 and 400 Myr); the grey intensities correspond to the logarithms of particle numbers in the pixels. Middle and 
bottom panels show the dependences of various disc parameters on the cylindrical radius R for various moments of time. Here n is the 
number of particles in concentric cylindrical layers; 1Z\ji is twice the median of the value \z\ (it is a parameter of the disc thickness, see 
ISotnikova fc Rodionovll2006l) ; v<p, ctr, (Tm, <t z are four moments of the velocity distribution. 



Another important feature of the FSLC.O model is 
that the velocity distribution in this model is far from the 
Schwarzschild distribution. The velocity distributions in the 
solar neighbourhood (near R = 8 kpc) are shown in Fig. [5] 
Initially, both radial and azimuthal velocity distributions 
are far from Gaussian, although such unusual distributions 
are more or less in equilibrium. At least, they are conserved 
during the initial stages of evolution. However, these dis- 
tributions are unstable, and they change substantially after 
about 500 Myr. After about 1 Gyr, the distributions are 
smoothed and tend to Gaussians. 

Above we have discussed the self-consistent evolution 
of the FSLC.O model. It is also interesting, however, to ex- 
amine a non-self-consistent evolution of this model; that is, 
to calculate the evolution of N test particles in the total 



potential of the FSLC model. As expected, the model prac- 
tically does not change during such "evolution" (at least on 
a time-scale of 10 Gyr). In particular, the density profile 
and non-Schwarzschild velocity distribution do not change. 
Thus, this shows again that this non-Schwarzschid velocity 
distribution is in equilibrium. 

All models constructed using the Orbit. NB approach 
have the following properties. They are close to equilibrium. 
Moreover, the velocity distributions in the models are non- 
Schwarzschild ones and may have various forms. However, 
the velocity distributions tend to the Schwarzschild distri- 
bution during the dynamical evolution on a time-scale of 
1 Gyr. Thus, although these models are close to equilib- 
rium, they are probably non-physical because of their non- 
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Figure 4. Comparison of profiles of the velocity distribution moments calculated from the Jeans equations and from the FSLC.O model. 
All moments for the disc were calculated inside the region \z\ < 0.1 kpc. Left panel: the solid line corresponds to the value v v for the 
model, and the dashed line corresponds to the same value calculated from the Jeans equation (the first equation of the system J9}), 
where the values an and a v were taken from the model. Middle panel: the solid line corresponds to the value a v for the model, and 
the dashed line corresponds to the same value calculated from the Jeans equation (the second equation of the system (0), where the 
values v v and an were taken from the model. Right panel: the solid line corresponds to the value a z for the model, and the dashed line 
corresponds to the same value calculated from the Jeans equation (the third equation of the system 10). 



Schwarzschild velocity distributions. The arguments are as 
follows. 

• The velocity distribution of the stars in the solar neigh- 
bourhood is similar to the Schwarzsch ild distribution (see, 
for example. [Binnev fc Merrineldlll99ct ). 

• The constructed non-Schwarzschild velocity distribu- 
tions are almost in equilibrium, but are unstable. The ve- 
locity distributions always tend to the Schwarzschild one 
during the evolution. 

• In the models constructed using the Nbody.NB ap- 
proach, the final velocity distributions are close to the 
Schwarzschild distribution (see Fig. [6}. We note that the 
Nbody.NB and Orbit. NB approaches differ only in the 
method of evolution simulations in the iterations (see Sec- 
tion 

Generally speaking, we can assume that such unusual 
non-Schwarzschild velocity distributions will survive only in 
the "hothouse" conditions of the Orbit. NB approach be- 
cause the evolution simulation in this approach is carried 
out non-self-consistently (see Section T3.2,3p . When the con- 
ditions are close to realistic (as in the Nbody.NB approach), 
such distributions are smoothed and gradually converge to 
the Schwarzschild distribution. 



At the same time, the models constructed by the 
Nbody.NB approach are the same (within the noise level) 
for arbitrary initial states at the same pdisk(R, z), $ ex t(i?, z) 
and fixed L z . Moreover, the velocity distributions in con- 
structed models are always close to the Schwarzschild dis- 
tribution. Based on this fact and the probable non-physical 
character of the Orbit. NB models, one can formulate a hy- 
pothesis that the "physical" discs in the equilibrium state 
are unique. This hypothesis could be formulated as follows. 
When the functions Pdisk(^, z) and "3> ex t(fi, z) are fixed and 
a fraction of the kinetic energy in the residual motions is 
also fixed (e.g. the value of the angular momentum L z is 
fixed), not more than one physical model in the equilibrium 
state may exist (the physical model is the model that can 
exist in reality). We think that the main feature of a real 
stellar disc is an almost Schwarzschild velocity distribution. 

Whether or not our hypothesis is true is very impor- 
tant. If we know the mass distribution in the Galaxy and 
our hypothesis is valid, we can use our Nbody.NB method 
to reconstruct the velocity distribution in the Galactic disc. 
In other words, we can construct a phase model of the Galac- 
tic disc, and this model will correspond exactly to the real 
Galactic disc. 



4.2 Uniqueness hypothesis 

In RS06, the authors formulated a hypothesis of uniqueness: 
not more than one equilibrium model (one equilibrium dis- 
tribution function) can exist for a fixed density Pdisk(-R, z) 
and potential "J> cx t(i?, z) and a fixed kinetic energy fraction 
of residual motions (e.g. a fixed angular momentum L z ). 

We can now say that in such a form the hypoth- 
esis is false. We can construct using the Orbit. NB ap- 
proach as many equilibrium models as we want for the same 
Pdisk{R, z), $ext(-R, z) and fixed L z . However, although these 
models are close to equilibrium, they probably bear no rela- 
tion to actual stellar systems. 



5 PHASE MODELS OF THE GALACTIC 
STELLAR DISK 

5.1 Choice of model among the family of models 

For the construction of phase models of the Milky Way 
Galactic disc, we use the Nbody.NB approach (see Sec- 
tion 13. 2p . Using this method, one can construct a family 
of models for the fixed functions Pdisk(R,z) and &e*t(R, z). 
This is a one-parameter family, and the parameter is the 
fraction of kinetic energy of residual motions (in other words, 
the disc "heat" degree). In our approach, this parameter is 
the value of L z . We emphasize that, using the Nbody.NB 
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7.5 kpc < R < 8.5 kpc 




Figure 5. The velocity distributions for the FSLC.O model at various moments of time (0, 100, 200, 600, and 1000 Myr). The region 
of the disc within the range 7.5 kpc < R < 8.5 kpc is considered. Left column: two-dimensional velocity distribution (the abscissa is the 
radial velocity vr and the ordinate is the azimuthal velocity v v )\ the grey intensities correspond to the numbers of particles that have 
velocities in the corresponding pixels. Middle column: one-dimensional distribution of the velocity vr. Right column: one-dimensional 
distribution of the velocity v v . In the middle and right columns, the solid lines show the model distributions, and the dashed lines 
correspond to the Gaussians with parameters (mean and dispersion) taken from the models. 
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Figure 6. The distributions of radial and azimuthal velocities in the solar neighbourhood in the models constructed using the Nbody.NB 
approach. We consider the disc region of 7.5 kpc < R < 8.5 kpc. Left panels: the distributions of radial velocities. Right panels: the 
distributions of azimuthal velocities. The upper pictures were constructed for the FSLC.N model, and the lower ones for the DB2.N 
model. Description of both models is given in Section 15.21 The solid lines show the model distributions, and the dashed lines correspond 
to the Gaussians with parameters (mean and dispersion) taken from the models. 



approach, when we fix the value of L z we obtain the same 
model (within the noise level) regardless of the initial state. 

A family of the JV-body models for the FSLC density 
model is shown in Fig. [7J as an example. The parameters 
of the models are given below. One can ask how we choose 
the best-fitting model among the family. We have compar- 
atively reliable kinematic Galact ic parameters in the solar 
neighbourhood (see, for e xample, iBinnev fc~ Mcrrificldl ll998l : 
iDehnen fc Binney|[l998bT ) , so it is reasonable to use them to 
choose the model. One needs to choose any one parameter 
that on the one hand is well known in the solar neighbour- 
hood, and on the other depends strongly on the disc "heat" . 
In other words, this value has to be strongly dependent on 
the value of L z . For example, the velocity ellipsoid parame- 
ters v v , or, Otpi and a z are well known. The value a z is not 



suitable, because it does not depend on L z (see Fig. [7] and 
RS06). The value v v is also unsuitable, because it depends 
weakly on L z for cold models (see Fig. [JJ. The choice be- 
tween the two values or and <r v is somewhat arbitrary. We 
prefer the value or for the model choice. 

There are various estimates of the value of or in the 
literature. We have chosen the value or — 35 kms -1 that 
was estimate d using an extrapo lation to the zero heliocentric 
distance (see lOrlov et al.l 2006). In addition, we adopted the 
solar distanc e from the Galactic centre as R® — 8 kpc (see, 
for example, iNikiforovlliooi ; lAvedisovalkoOSl ). 

As a result, we have chosen the model in which the ra- 
dial velocity dispersion in the solar neighbourhood is about 
35 kms -1 . As the solar neighbourhood, we have chosen the 
region 7.9 kpc < R < 8.1 kpc and \z\ < 0.1 kpc. 
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Figure 7. The family of phase models constructed using the Nbody.NB approach for the FSLC density model. For each model, we 
show the dependences of v v , an, (Tip, and a z on R. The moments were calculated in cylindrical layers infinite along the z-axis. In the 
upper left panel, the circular velocity curve is shown by the thick solid line. Models constructed for L z = 0.302, 0.3, 0.29, 0.28, 0.27 are 
represented. The model with L z = 0.302 is the FSLC.N model corresponding to the Galactic stellar disc. 



5.2 Models FSLC.N and DB2.N 

We have considered two families of stellar disc models con- 
structed for two Galactic density models (FSLC and DB2). 
The families of phase models were constructed using the 
Nbody.NB approach. From every family, we have chosen the 
stellar disc model that corresponds to the Galactic disc in 
the solar neighbourhood (in terms of the radial velocity dis- 
persion). Hereafter we refer to these models as FSLC.N and 
DB2.N. 

The family of models for the FSLC density model is 
shown in Fig. [7] It was constructed in the following way. 
We took as pdisk the disc density distribution from the 
FSLC model, and as $ cx t the potential arising from all 
FSLC model components except the disc. The disc den- 
sity was taken as equal to zero at R > i? m ax = 30 kpc 
or \z\ > z max = 10 kpc. We considered an initial cold model 
in which the particle orbits are circular. Four iteration sets 
with increasing accuracies were calculated. The parameters 
of the sets are shown in Table [3] The number of neighbours 
in the velocity distribution transfer is n n b = 10. The time of 
one iteration is U = 50 Myr. The FSLC.N model was con- 
structed for the angular momentum L z = 0.302 (this value 
is given in our system of units as described in Section [4. 1|) . 

The family of models for the DB2 density model, in par- 



Table 3. The parameters of four iteration sets for the family of 
models constructed using the Nbody.NB approach for the FSLC 
and DB2 density models. Here, n- lt is the number of iterations, N 
is the number of bodies, dt is the integration step, ej is the soft- 
ening length for the FSLC model, and is the softening length 
for the DB2 model. The softening length was ch osen using the 
recommendations of lRodionov &: Sotn ikova ( 20051). 
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ticular the DB2.N model, was constructed in the following 
way. The density of the thin stellar disc in the DB2 model 
was taken as pdisk, and the potential arising from all com- 
ponents except the thin stellar disc was taken as <& ox t. The 
disc density at R > i? max = 30 kpc or \z\ > 2 ma x = 10 kpc 
was taken as equal to zero. This is the same condition as 
in the FSLC model. Here we again consider the cold initial 
model with circular orbits. Four sets of iterations were again 
calculated. The parameters of these sets are shown in Table 
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3. We have adopted the parameters n n b = 10, U = 50 Myr. 
The DB2.N model was constructed for the angular momen- 
tum L z — 0.1595 (this value is given in our system of units 
as described in Section [4TTJ) . 

The radial profiles of the velocity distribution moments 
for the FSLC.N and DB2.N models are shown in Fig. 
Let us consider the profiles of an and a v . In the central 
parts of the models, these profiles are very different. This 
is caused by the more massive and concentrated bulge in 
the FSLC.N model with respect to the DB2.N model (see 
Fig. [T] in Section [2]|. However, the profiles of an and a v for 
the two models are similar from about 2 — 3 kpc. This is 
in spite of the difference between the initial density models. 
In particular, one can observe a different relative disc con- 
tribution in the total mass and rotation curve (see Fig. [TJ. 
Fig. [5] also shows that the profiles of a z are different. How- 
ever, we note that the value of a z at any point is defined 
only by the density distribution. This fact is a consequence 
of the last Jeans equation . Our phase models satisfy the 
Jeans equations very well, however, and therefore the differ- 
ences in the a z profiles are explained by the differences in 
the density models. 

The initial stages of the evolution for the FSLC.N and 
DB2.N models are shown in Figsl9land llOI Both models con- 
serve the structural and dynamical parameters in the early 
stages of the evolution very well. Therefore both models are 
close to equilibrium. 



6 CONCLUSIONS 

In this paper, we have discussed the problem of construct- 
ing a phase model of the Galactic stellar disc. We used 
the iterative method proposed in RS06. The realization of 
the iterative method described in RS06 (Nbody.SCH ap- 
proach) has some disadvantages. The Nbody.SCH method 
has a problem with the construction of a relatively cold 
stellar disc. Furthermore, the Nbody.SCH method cannot 
be directly applied to stellar systems with arbitrary geom- 
etry. The main goal of this study was to develop a new 
realization of the iterative method without the disadvan- 
tages of the Nbody.SCH approach. For this purpose, we con- 
sidered a number of modifications of the iterative method 
(Nbody.SCH, Orbit.SCH, Nbody.NB, and Orbit.NB). The 
Nbody.NB method satisfied our conditions. This method can 
be directly applied for the construction of phase models of 
stellar systems with arbitrary geometry. The method is sim- 
ple in terms both of understanding and of implementation. 
Using the Nbody.NB approach, we have constructed phase 
models of the Galactic disc for two realistic density models 
(suggested by | Flynn, Sommer-Larsen fc ChristensenI Il99rj ; 



1M 

iDehnen fc Binnevlll998ah . 

For a given mass distribution model of the Galaxy, we 
can construct a one-parameter family of equilibrium models 
of the Galactic disc. From this family we can choose a unique 
model using local kinematic parameters that are known from 
the Hipparcos data. There are, however, two important ques- 
tions. Is there an equilibrium disc model besides the models 
from this one-parameter family? Does the real Galactic disc 
belong to this one-parameter family? (Of course this family 
should be constructed for a real mass model of the Galaxy.) 
The answer to the first question is affirmative. Using the Or- 



bit.NB method, we can construct many models besides the 
models from this one-parameter family. We have, however, 
shown that these models are probably non-physical. Based 
on this fact, we suppose that all models except models from 
this one-parameter family are non-physical. Consequently, 
we think that the answer to the second question is also af- 
firmative. However, we cannot strictly prove it as yet. We 
think that the key to the proof is Schwarzschild velocity dis- 
tributions. The models from our one-parameter family have 
almost such velocity distributions, and we think that the 
velocity distributions in real galactic discs are close to the 
Schwarzschild distribution. 

We now discuss possible applications of our Nbody.NB 
method. 

This method can be used to construct phase models of 
stellar systems with arbitrary geometry. For example, the 
method can be used to construct phase models of triaxial 
elliptical galaxies. 

Here we have used the iterative method to construct 
a phase model of the stellar disc of a spiral galaxy. In the 
future, we are going to construct a self-consistent model of a 
spiral galaxy (including live halo and bulge) using a modified 
Nbody.NB approach. 

Another, rather more important, direction of future in- 
vestigations is a comparison of our iterative models with 
observations of real galaxies. This comparison will make it 
possible to derive from observations the constraints on un- 
observable parameters of galaxies (for example dark matter 
mass and its profile). Of course, such a comparison should 
first be carried out for the Milky Way. One of the possi- 
ble schemes is as follows. If we know the mass distribution 
in the Galaxy then using our method we can construct a 
phase model of the Galactic disc. The mass distribution in 
the Galaxy contains two parts — visible and dark matter. 
From the phase model of the Galactic disc, we can derive the 
profiles of stellar kinematic parameters (profiles of velocity 
dispersions and mean azimuthal velocity). Let us assume 
that we know from observations the mass distribution of 
visible matter and the profiles of stellar kinematic param- 
eters. Then, using the iterative method, we can put some 
constraints on the dark matter mass distribution. The dis- 
tribution of the dark matter should be such that the iterative 
model has the observable kinematics. At the moment, obser- 
vational profiles of stellar kinematic parameters have fairly 
large uncertainties; however, it is expected that the GAIA 
mission will provide a much higher accuracy for these data. 
We are planning to study which constraints on the Milky 
Way parameters we can obtain from GAIA using our itera- 
tive models. 

The iterative method can also be used to derive the ve- 
locity dispersion profiles from observations. It is impossible 
to obtain the velocity dispersion profiles for external spi- 
ral galaxies directly from observations, which can yield only 
the line-of-sight velocity, V\ os , and the line-of-sight velocity 
dispersion, cti os . However, we can use the Jeans equations 
to derive the velocity dispersion profiles from the observed 
quantities v\ OE and <ti os . In practice, however, one has to in- 
clude some a priori assumptions about the form of the veloc- 
ity dispersion profiles. In the literature, authors have made 
slightly different assumptions, but all of them are based on 
the hypothesis that the velocity dispersion in the radial di- 



14 S.A. Rodionov, V.V. Orlov 




rection is proportional to the velocity dispersion in the ver- 
tical direction (see section 3 in RS06). 

Using our iterative method, we can derive the velocity 
dispersion profiles from line-of-sight parameters without any 
additional assumptions. The general idea is as follows. From 
observations we have (more or less precisely) the distribution 
of visible mass and some constraints on the distribution of 
dark matter. If we fix the mass distribution in a galaxy (for 
visible and dark matter) then, using the iterative method, 
we can construct a one-parameter family of phase models of 
the galactic disc. The parameter of this family is the degree 
of disc heating (for example the fraction of the kinetic energy 
of the disc contained in residual motions). Using observable 
profiles of Vioa and <Ti oa , we can choose a unique model from 
this family. As a result, we will have a phase model of the 
galactic disc. From this model we can calculate the velocity 
dispersion profiles. 
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Figure 9. Initial evolution stages of the FSLC.N model. The same values are shown as in Fig. \3\ The evolution simulation used the 
following parameters: number of bodies N = 100 000, integr ation step dt = 0.04 Myr, so ftening length e = 0.025 kpc. The two last 
parameters were chosen according to the recommendations of lRodionov fc Sotn ikova ( 2 0051) . 
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Figure 10. Initial evolution stages of the DB2.N model. The same values are shown as in Fig. \3\ The evolution simulation used the 
following parameters: number of bodies N = 100 000, integr ation step dt = 0.04 Myr, so ftening length e = 0.015 kpc. The two last 
parameters were chosen according to the recommendations of Rodionov & Sotnikova (2005). 



